% This is a test file that solves
% "minimize f(X)" with rank(X) = r
%
% Refer "Low-rank optimization with trace norm penalty"
% Technical report (arXiv 1112.2318), 2011 for further details.
%
% Authors:
% Bamdev Mishra <b.mishra@ulg.ac.be>


clear all; close all; clc;
RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));

cd ..

%% Setup
d1 = 5000;
d2 = 5000;
r = 5;
over_sampling = 5;
lambda = 0;

dof = (d1 + d2 - r) * r;
f = over_sampling * dof / (d1 * d2);
fprintf('Rank %i matrix of size %i times %i and OS = %i\n', r, d1, d2, over_sampling);

% Options
verbosity = true;
random_initialization = false;
compute_predictions = false; % Compute predictions at each iteration

%% Optimization algorithm is the trust-region algorithm
optim_algo = @tr_matrix_completion;



%% Paramters
params.tol = 1e-12; % Absolute tolerance for Trust-region algorithm
params.vtol = 1e-12; % Relative tolerance for Trust-region algorithm
params.grad_tol = 1e-5; % Absolute tolerance for Trust-region algorithm
params.dg_tol = 1e-5; % Absolute tolerance for duality gap
params.dg_vtol = 1e-5; % Relative tolerance for duality gap
params.smax_tol = 1e-3; % Absolute tolerance for duality gap
params.max_iter_tr = 50; % Maximum iterations for the TR algorithm
params.verb = verbosity; % Dislpaly output if asked
params.compute_predictions = compute_predictions; % Compute predictions each iteration if asked


params = default_params(params);


%% Generating random data
data = generate_randn_data(d1, d2, r, f);


% Testing data 
data.ts.nentries = data.ls.nentries;
data.ts.rows = randi(d1, data.ts.nentries, 1);
data.ts.cols = randi(d2, data.ts.nentries, 1);
data.ts.entries = partXY(data.Gs', data.Hs',data.ts.rows,data.ts.cols,data.ts.nentries)';


%% Creating a sparse structure
sparse_structure = sparse(data.ls.rows,data.ls.cols,data.ls.entries,d1,d2);
sparse_structure2 = sparse(data.ls.rows,data.ls.cols,2*data.ls.entries,d1,d2);



%% Setup initial model
model.d1 = d1;
model.d2 = d2;
model.p = r;
model.sparse_structure = sparse_structure; % sparse skeleton 1
model.sparse_structure2 = sparse_structure2; % sparse skeleton 2

model.lambda = lambda;



%% Functions
fun_set=@functions_nsym_polar_geometry; % File that specifies geometry of the search space UBV^T
fun_obj=@functions_matrix_completion_polar; % File that computes cost function, gradient of the matrix completion problem




%% Intialization
if ~random_initialization,
    [model.U model.B model.V] = svds(sparse_structure, r);
    fprintf('**** Initialization by taking %i dominant SVD\n', r);    
else
    G = randn(model.d1, model.p); H = randn(model.d2, model.p);
    [Qg Rg] = qr(G, 0);
    [Qh Rh] = qr(H, 0);
    [q1 b q2] = svd(Rg * Rh');
    model.U = Qg * q1;
    model.V = Qh * q2;
    model.B = b;
    fprintf('**** Random initialization\n');
end



%% Trust-region algorithm for fixed rank

tic;
[model_ubv, infos_ubv] = tr_matrix_completion(data.ts, data.ls, model, params);
toc;






